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Abstract 

We study the phase transitions of a restricted solid-on-solid model coupled 

to an Ising model, which can be derived from the coupled XY-Ising model. 

There are two kinds of phase transition lines. One is a Ising transition line and 

the other is surface roughening transition line. The latter is a KT transition 

line from the viewpoint of the XY model. Using a micro canonical Monte 

Carlo technique, we obtain a very accurate two dimensional phase diagram. 

The two transition lines are separate in all the parameter space we study. 

This result is strong evidence that the fully frustrated XY model orders by 

two separate transitions and that roughening and reconstruction transitions 

of crystal surfaces occur separately. 
PACS numbers:64.60.-i, 64.60. C, 75.40.Mg 
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I. INTRODUCTION 

Since the pioneering work of Teitel and JayaprakashS, two dimensional (2D) fully frus- 
trated XY (FFXY) models have been intensively studied over the last decade, both ana- 
lytically and numerically.lH! This model has a continuous U(l) symmetry associated with 
global spin rotations and also a discrete Z{2) symmetry because of the double degeneracy of 
the ground state vortex configurations. At sufficiently low temperature (T) the system will 
be completely ordered with true long range order in the Ising order parameter characterizing 
the vortex lattice and algebraic decay of correlations of the XY order parameter. As T is 
increased, the system will disorder by one of two possible scenarios. One is the simultaneous 
loss of both XY and Ising order and the other is a two stage process in which XY order 
is destroyed at a lower temperature than the Ising order. Since the original paper on this 
model,0 an enormous effort has gone into trying to clarify the nature of the phase transi- 
tion^) of this model and several related models believed to be in the same universality class 
have been studied numerically.lH! However, the nature of the transition(s) in the FFXY 
model and its generalizations is still somewhat controversial. Monte Carlo (MC) simula- 
tions of the Coulomb gas with half integer chargesi@ generally show a vortex lattice melting 
transition and a KT-like transition at slightly different temperatures which is consistent with 
the second scenario. However, the melting transition is expected to be in the Ising univer- 
sality class but the simulations yield exponents which are inconsistent with the Ising values. 
Moreover, at the KT-like transition the discontinuity in the helicity modulus is much larger 
than the universal value of 2/7ri. Yet another model is the coupled XF-Ising model! which 
is a coarse-grained version of the FFXY model. It was suggested that this model displayed 
both possible scenarios depending on the parameters of the model. Recent studies!! of the 
FFXY model with a nearest neighbor cosine interaction show that the disordering is by 
the second scenario with a KT transition followed by an Ising transition at a slightly higher 
temperature. From a theoretical point of view, the FFXY model may have separate XY 
and Ising transitions if the critical exponents of the vortex lattice melting are consistent 
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with the 2D Ising values and the jump in the helicity modulus at the XY transition has the 
universal value of 2/ir but most numerical studies are not consistent with these. However, 
in his recent numerical study, OlssoJ pointed out that, for conventional finite size scaling 
to hold, it is essential that the system size L is much larger than any other length scale. 
Since the vortex lattice melting temperature T& is only slightly larger than Tkt, the relevant 
length is £+(T c i), the XY correlation length at the lattice melting temperature. As is well 
known, this remains rather large quite far above Tkt and, by taking this into account, a 
good fit to the double transition scenario was obtained^ with a conventional KT transition 
followed by an Ising transition. 

Whatever the precise nature of the transitions in any model, an essential first step is to 
determine whether the system has two separate transitions or a single transition. It is of 
some interest from the point of view of conformal field theories to find a non-trivial model in 
which the central charge c > 1. Recently, a related class of models for transitions on crystal 
surfaces with competing reconstruction (Ising) and roughening (XY) degrees of freedom 
has been studied0@. Den Nijs has studied competition between surface roughening and 
reconstruction and tentatively concluded that the roughening and deconstruction transition 
lines meet and finally merge together into a single roughening induced simultaneous recon- 
struction transition line.0 On the other hand, Carlon in a recent study of the staggered 
six vertex modeljll has suggested that the two transition lines never merge but become 
very close to each other which is strong evidence for two separate transition in the FFXY 
model. In order to clarify whether the Ising and XY transition lines merge we introduce 
the restricted solid-on-solid (RSOS) model coupled to the Ising model. This RSOS-Ising 
model can be derived from the XY-lsing model by making a duality transformation on the 
XY degrees of freedom which introduces integer height variables on the dual lattice. The 
resulting RSOS-Ising model has both Ising and height (XY) degrees of freedom and, for 
this reason, can be considered to have the same features as the staggered six vertex model 
and reconstructed surface models. 

We investigate the RSOS-Ising system by a rejection free microcanonical MC technique^ 



which is appropriate to this system and also highly accurate. We have obtained an accurate 
phase diagram in a two parameter space and some critical exponents. As a result we have 
found there are Ising and roughening critical lines which never merge in all parameter space. 
This result is strong evidence against the single transition scenario in the fully frustrated 
XY model. 

In the following section we give a derivation of the RSOS-Ising model from the coupled 
XF-Ising system. In section III we present a brief explanation of the MC procedure^ as 
adapted for this system, methods of data analysis and MC results. Finally we present our 
conclusions. 



II. RESTRICTED SOLID-ON-SOLID MODEL COUPLED TO ISING MODEL 

From renormalization group and universality arguments the isotropic FFXY model can 
be described by a coupled XY-lsing model of the form 

- (3H = A J2( 1 + Wj) cos ( d i -6j) + c oYl a i a h (!) 

(hi) 

where 0, is the phase of the continuous spin at zth site and Oi = ±1 describes the chirality. 
The physical range of the phenomenological parameters is Aq > and Aq + Co > 0. This 
model has been intensively studied by Lee et aM By studying the behavior of the chiral 
(Ising) degrees of freedom <7j, they found a line of continuous transitions with a segment 
of non-Ising like critical exponents, in particular uj « 0.8 < 1 and r\i > 1/4. Based on 
this, they argued that the Ising and XY critical lines must have merged into a single line 
of simultaneous ordering of both degrees of freedom. However, they made no explicit study 
of the XY degree of freedom for technical reasons, which weakens their argument. To 
simultaneously study the XY and Ising critical lines, we make a duality transformation^ 
on the XY variables of Eq.(|l|) to obtain a restricted solid-on-solid model coupled to an Ising 
model, in which only discrete variables are involved, as described below. 
The partition function of the Xy-Ising model of Eq.(|l]) reads 
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where {a} denotes the sum over all configurations of Ising spins, N is the number of lattice 
sites and Kij = Aq(1 + <Ji<jj). Using the periodicity of the integrand in 9ij = 8{ — 9j we have 



exp [K i:j cos(9i - Bj)\ = 51 ex P [ in »i%] ex P ^( 



(3) 



7lij = — OO 



where exp V^n^) is a Fourier component of exp [K^ cos(0j — This yields 
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The integrations over the can be easily done and, noting that n»j = — riji, give iV con- 
straints on the 2N integers so the partition function has the form 
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where the A(z) denote the constraints at site i. For example, the shown in Fig. 1 

gives A(0) = 8n 01 +no2,n 3 o+n4 - The N constraints are automatically satisfied by introducing 
the integer height variables hj on the dual lattice shown in Fig. 1 as 

n i = h 1 - h 2 , n 02 = h 2 - h 3 
n 3 o = h 4 - h 3 , n 40 = h x - h 4 . 

In terms of the heights hi the partition function has the apparently simple form 
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where I and J are dual lattice sites and —oo<hj<oo. To this point everything is exact but 
Eq.(||) is not suitable for simulation because of the infinite range of integer height differences 
hj — hj. The Fourier component is 

1 r 27T 



exp 



V(hj - h d 



2tt Jo 



exp [—i(hi - hj)9ij] exp [A (l + er^-) cos 9^) d%, (7) 
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where the nearest neighbor bonds 13 on the dual lattice and ij on the original lattice cross 
each other. If a^aj = — 1, the Fourier component exp V(h) = 8^. A step in hi cannot 
cross an Ising domain wall. When a^aj = 1 Eq. (|7]) becomes Ih(2Ao), where Ih is a Bessel 
function of order h with imaginary argument and h = hj — hj and //i(0) = 5h,o- From 
the asymptotic behavior of the Bessel functions we have Ih(2A )/ I (2A ) « exp(— h 2 /4:A ), 
which decays rapidly as h increases if the inverse temperature A Q is not too large. Therefore, 
it is a reasonable approximation to keep only the contributions from hi — hj = 0, ±1 when 
OiOj = 1. 

In view of these arguments, an effective Hamiltonian for the system is 

-/3H = CY,<n<rj-A y Z \hi~hj\ (8) 

where (hj — hj) = 0, ±1 if the corresponding o^Oj = 1 and hi = hj otherwise. By comparing 
the Boltzmann weights derived from the Hamiltonians (|I|) and (^j) we obtain the relations 
exp(-A) = h(2A ) / I (2A ) and exp(2C) = exp(2C )/o(2A ). A is inversely related to A 
and C is a linear sum of Cq and A . Eq. (§) is the Hamiltonian for a RSOS-Ising model and 
the heights hi and Ising spins <7j are coupled via the constraint that a step in h is forbidden 
to cross a domain wall in a. 



III. MONTE CARLO ANALYSIS 

Due to the restrictions on the spin and height variables discussed in the last section, some 
spin-height configurations are not allowed which we call forbidden configurations. In our 
MC procedure, we select a spin or height variable randomly and update the configuration 
by a trial move. The trial move for the spin variable a is a simple flipping, i.e., a — > —a. 
For the height variable h, the trial move is either raise or lower the selected height by 1 unit, 
i.e., h — > h ± 1. If the attempted move takes the system to a forbidden configuration, we 
call the move forbidden and otherwise we call the move allowed. The number of forbidden 
moves changes according to the configuration of the neighboring spins and heights. For 



6 



accuracy and efficiency, we use the microcanonical MC methodlll introduced by one of the 
authors to study thermodynamic functions of the system in the 2D parameter space (A, C). 
We further apply a rejection-free technique to speed up the MC procedure.^ To apply 
the microcanonical MC method to the RSOS-Ising model we first introduce an energy-like 
quantity 

E = E l + E R , (9) 

where E\ = — J2a,j) Vi&j an d -Er = 2J2n,j) \hi — hj\. The factor 2 in Er is introduced for 
convenience to make the minimum energy changes AE-& = AEj = 4 for both E-& and Ej. The 
quantities to be computed are the number of configurations u(E, E{) and the microcanonical 
average Q(E, E{) of a thermodynamic quantity Q at a set of values of E and Ej. 

Since the rejection free MC technique^! is not well known, we present a brief account of 
the microcanonical MC method as adapted to the RSOS-Ising system. We set up a random 
walk through a configuration space restricted to a small energy square in the (E, E{) domain 
given by 

E® - qAE/2 <E< E {i+1) + qAE/2 (10) 
E[ j) - qAEi/2 <Ej< E[ i+1) + qAEi/2 (11) 

where AE = 4 is the minimum energy separation in E, q is the coordination number of the 
lattice, E^ = E® + AE and E[ j+X) = E[ j) + AE h qAE/2 and qAE l /2 are the maximum 
energy changes in E and E\ when a spin or height variable is updated to a new state. The 
size of the square given by Eqs. (|10D and ([11]) is chosen to make any allowed move keep 
the new configuration within this square if the current configuration is in the innermost 
square, a configuration whose E and E\ values belonging to one of four sets, (E^,E^ j) ), 
(£7( < + 1 ),Ep ) ), (E^\Ei 3+1) ), and (E^ i+1 \ E^ +1) ). 

Since there are N candidates for the trial spin move and 2N candidates for the trial 
height move, there is a total of iV = 3iV of candidates for the trial move. Then out of 
this N, there are iV a i wd allowed and iVfbdn = N — iV alwd forbidden moves. We tabulate both 



allowed and forbidden moves in the same look-up table. The allowed moves are divided into 
two groups, one consisting of acceptable moves which keep the walker within the predefined 
energy domain given by Eqs. ( PD and ([II]) and the other of rejectable moves which take 
the walker outside this square. If the walker is initially at a configuration belonging to the 
innermost square, all allowed moves are acceptable. On the other hand, if the walker is not 
at configuration belonging to the innermost four energy values, then some allowed moves are 
acceptable and some are rejectable since some moves will take the random walker outside 
the square of Eqs. ( |lOl) and (|TT|) as the largest energy changes in E and Ei are qAE/2 and 
qAEi/2. 

Whenever the walker visits one of the configuration belonging to the innermost square we 
count the number of visit and take samples of Q's. Next we select a trial move randomly from 
all possible N moves. If the selected move is an allowed one, we update the configuration 
since it is always acceptable. If the selected move is forbidden we discard the move but count 
the current configuration once more and sample again for the same configuration to satisfy 
the condition of detailed balance. This procedure is equivalent to assigning the correct 
weight to the sampling configurations. When the walker is at a configuration not belonging 
to the innermost four energy values, we make a random selection of trial move from the group 
of acceptable moves only and let the walker move to a new configuration without wasting 
computing time. This is possible since we are not taking data at this configuration. In 
this way we can estimate the number of visits for the innermost four energy configurations, 
which give relative ratios E^ +1) )/uj(E^ , E[ j) ) and u(E^ i+1 \ E^)/uj(E^, E[ j) ) for all 

i and j. The sampled data for Q will furnish microcanonical averages Q(E,E{). Since the 
exact number of configurations u(E, E{) for small values of E and E\ are easily computed, 
we can recursively estimate all tv(E, Ei)'s from these. To obtain more accurate data we use 
a two-step MC simulation. First we accurately estimate uj(E^ l+lS) ) /u(E^) for all i, where 
Q(E^) = We then estimate Q(E^\E^ +1) )/u(E^\eP) for all j. Now we 

can easily estimate all u(E, E{) consistently. The microcanonical average is 
also estimated directly from data taken during this MC process. To check the validity of our 
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procedure we compare the estimated uj(E r , 0) and u(Ei, Ei) for all system sizes with those 
of the pure RSOS model and those of the pure Ising model respectively. We calculated the 
number of configurations for the pure RSOS model independently and the those of the pure 
Ising model are well knownS Both of them agree well with each other to within expected 
statistical deviations. 

The canonical average of any thermodynamic quantity Q is obtained by 

(Q) = Eee^m, Q(E, E h M R , M{)u(E, E h M R , MQ expj-CEj - AE R /2) 

W; T.ee^mME^i^^^^M-CEj-AE^) ' 1 ' 

where Mi = L~ 2 \J2i&i\ is the Ising magnetization, Mr is a RSOS magnetization to be 

defined in Eq.fJfJ]) below and Q{E, Ej, Mr, Mi) and u(E, Ei, Mr, Mi) are the microcanon- 

ical values of the thermodynamic quantity Q and the number of configurations for given 

(E, Ei, Mr, Mi), respectively. In Eq.(|12|), the exponent of the Boltzmann weight factor does 

not involve any magnetization. We thus estimate only uj(E, Ei) = J2m r Mi ^(E, Ei, Mr, Mi) 

and Q(E, E{) = £ MrMi u(E, E h Mr, M x )Q{E, E h Mr, M^/Cj^E, El). If Q depends only on 

E and E\, as in the case of the specific heat, we need only u(E, E{) since Q(E, E\) = 

Q(E,Ei). Otherwise we need to calculate both uj(E,Ei) and Q(E,Ei). The canonical 

average of Q can be written as 

, n) Zee, Q(E, Ej)u(E, Ej) exp(-CE l - AE R /2) 

W; Ee^u&Edexpi-CEk-AEzffl ' 1 } 

All thermodynamic quantities we use in this paper depend only on the energies E and Ej 
or only on magnetizations Mr and Mj. For example, Q(E, Ex, Mr, Mi) = E for the total 
energy, Q = Mi for the Ising magnetization and Q = Sm^m for the magnetization probability 
distribution P^(M). In all cases, samplings are performed in the 2D plane of (E,E{). The 
number of possible E and E\ is proportional to L 2 and the simulation time increases as 
L A where L is the lattice size. This makes it difficult to simulate systems of large size. 
Nevertheless, once Q(E,Ej) and uj(E,Ei) are estimated, the canonical average (Q) can be 
easily calculated for all A and C from Eq. ([13]). 

As the specific heat of the 2D XY model has no diverging peak as L — > oo, the RSOS 
part of the specific heat may also not diverge which implies the transition temperature 
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may not be easy to determine accurately. For the 2D XY model, which has vanishing 
magnetization in the bulk limit for any finite temperature, the helicity modulus, or stiffness 
constant which is a measure of the quasi long-range order of the order parameter describing 
the continuous 17(1) symmetry and the superfhud density of a Bose superfluidB, may be used 
to determine the XY transition temperature from the finite size scaling formS at the KT 
point. The analogous quantity in the RSOS model is the free energy of a step in the height 
which has the finite size scaling form@ F S (L) = ir/A + l/(A + BlnL) at the roughening 
temperature. However, the computation of this requires the partition function with two 
different boundary conditions which is time consuming. As an alternative method, we know 
that the roughening transition in the RSOS model is one between flat and rough phases. 
We introduce an order parameter, analogous to the magnetization of an Ising model, Mr, 
which measures the flatness of the RSOS system, as 

M R = L- 2 Y,™ R (h I ), (14) 
i 

where m R (hj) = 1 if the height hi is equal to the majority height and m R = otherwise. 
It is clear from this definition that Mr = 1 for a completely flat phase and vanishes in a 
rough phase in the thermodynamic limit. In Figs. 2(a) and 2(b) we plot Mr and M\ in 
the parameter space (A, C) for lattice size L = 22. It is easily seen from the figures that 
there are three different phases depending on the parameters A and C: Ising ordered rough 
(IOR), Ising ordered flat (IOF) and Ising disordered flat (DOF). The IOR and DOF phases 
are separated by the IOF phase in the low temperature region though it is not clear whether, 
in the region of small A and C, the two phases are also separated or not. The corresponding 
phases of the Xy-Ising model of Eq.fll]) are XY and Ising order (IOR), XY disorder and 
Ising order (IOF) and XY and Ising disorder (DOF). It is interesting to note that the fourth 
phase of XY order and Ising disorder, corresponding to an Ising disordered rough phase of 
the RSOS-Ising model, does not exist. 

We now examine the behavior of the specific heat. For simplicity, we consider Cr and 
Cj corresponding to the RSOS and Ising specific heats, respectively. These are defined as 
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L 2 C S = -D 2 d(E s ) /dD, where D,S= A,R and C,l respectively. In Figs. 3(a) and 3(b) we 
plot the RSOS and Ising specific heats for lattice size L = 22. We have obtained the specific 
heats for several lattice sizes up to L = 22. The RSOS specific heat ridge parallel to the 
C-axis ceases to grow with increasing L for L > 18 as for the standard XY model while the 
Ising specific heat ridge continues to grow as the system size increases. From the specific 
heats of several lattice sizes we find that the growth rate of the specific heat ridges in the 
small A and C region is no higher than that of the Ising ridge parallel to A-axis which implies 
that the exponent a < for either contribution. Because of the slow increase of the Ising 
ridge with L and the limitation to fairly small L, we are unable to obtain any reasonable 
estimate of the exponent a/u. 

Having established the existence of roughening and Ising transitions we can now locate 
the precise location of the two transition lines. It is relatively simple to locate the Ising tran- 
sition line by using the probability distribution of the Ising magnetization Mi, P L (Mi, A, C). 
In a phase with Ising order, Pl(Mj = 0, A, C) decreases to zero as the system size increases 
while Pl(Mj = 0, A, C) increases as L increases in an Ising disordered phase. We then ex- 
pect that there will be a unique crossing point (A*,C*) in a plot of Pl(M = 0, A, C) curves 
for several lattice sizes along a path in the A-C plane if the path crosses the transition line. 
Note that this does not use the concept of finite size scaling. This probability distribution 
method is a powerful tool to locate the transition point if the system is appropriate to this 
type of analysis. We have obtained Pl(Mi = 0, A, C) curves in the A-C plane for lattice sizes 
L = 10, 14, 18, 22 and found that there exists a well defined line of intersection points from 
several Pl(Mj = 0,A, C) surfaces which is shown in Fig. 4(b). The estimated positional 
uncertainty of the Ising transition line in the A-C plane is at most 0.001. Unfortunately, 
it is difficult to find any RSOS order parameter suitable for this type of analysis to locate 
the transition line and, as an alternative method, we use the position of the peak in the 
susceptibility of the RSOS order parameter M R . In the rough phase, this susceptibility is 
infinite in the thermodynamic limit and, for finite L, the susceptibility peak grows as L 2 ~ n . 
In the thermodynamic limit, the position of this diverging peak gives the transition point. 
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The peak position in a finite system, for fixed C, approaches its bulk value A* as 

\A peak (L) - A*\ ocL^, (15) 

where A peak (L) is the peak position for lattice size L. In Fig. 4(a) we plot A peak (L) versus 
L^ 1 for several C = 0.10, 0.15, 0.20, and 0.25 from bottom to top. All points are on almost 
straight lines which implies that Eq. (|1|) describes the scaling relation of the peak position 
with v — 1. Extrapolating to L^ 1 = 0, we obtain the bulk value of A* for several fixed 
C's. The specific heat exponent a < along the roughening transition line, as we argued 
earlier, and this implies v > 1. In fact, at the roughening transition point which is the KT 
point in the XY model, the exponent v = oo. If v > 1, A pe&k is convex up in L~ l and the 
extrapolated A* under the assumption v = 1 is an upper bound for the transition point. 
In Fig. 4(b) we draw the roughening transition line extrapolated from Eq. (|15D with v = 1 
together with the Ising transition line. The IOR and DOF phases are completely separated 
by the IOF phase in all parameter space (A, C). Since the phases are separated when we 
use an upper bound for A*, they will certainly be separate for the true values of A* as this 
will increase the separation. The phase trajectory of the FFXY model goes from IOR to 
DOF as the temperature increases from zero and the result is strong evidence for the double 
transition scenario in the FFXy model. 

Since the two phase transition lines are accurately located some critical parameters along 
these lines can be estimated. We have calculated the Ising susceptibility from the usual 
definition and estimated the exponent j/u along the Ising transition line using sizes L = 
6, 10, 14, 18, 22 for the fit (See Fig. 5). The exponent 7/V is consistent with the Ising value 
down to A = 0.8 but deviates significantly for smaller A where the two transition lines 
become close. Presumably this discrepancy is due to large corrections to scaling due to large 
fluctuations of the RSOS order and the limited system sizes. The height-height correlation 
function of the RSOS model defined by G(ti — rj) = ((hi — hj) 2 ) behaves as G(rj — rj) ~ 
a(A, C) log(|r/ — rj|) in the rough phase. The coefficient a(A, C) changes continuously in the 
rough phase and has a universal value 2/tt 2 at the roughening transition point .0 The height 
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width a 2 is defined by J2i(hi — h) 2 /L 2 , where h is the average height and, by integrating 
the height-height correlation function, a 2 as | log(L) + b in the rough phase. From the MC 
simulations, we have obtained a 2 (A, C) in the IOR phase and found that the height width a 2 
agrees well with the log L scaling form in the whole IOR region. The coefficient a increases as 
A decreases from A* into the rough phase, as expected. Along the roughening transition line 
where two transitions are not close to each other, the coefficient a ps 0.16, which is somewhat 
smaller than the universal value 2/n 2 ps 0.203. We find that the coefficient a 0.19 along 
the roughening transition line where two transitions are very close to each other. We do 
not use any scaling or extrapolation method to locate the RSOS transition points in the 
region where two transitions are well separated but simply use the RSOS susceptibility peak 
positions of the largest lattice size, L = 22, as the transition points. This peak position 
moves very slowly toward smaller A as the system size becomes larger, as expected. In fact 
the estimated a along the roughening transition line in that region decreases slightly as the 
system size increases. This means that the estimated transition line actually lies in the flat 
region close to the real roughening transition line which will lie to the left of our estimated 
line. We have located a more reasonable roughening transition line in that region using the 
criterion a = 2/ir 2 and find that the transition line shifts to the left about 6% in A. This 
is consistent with our arguments that the true A* is smaller than the value estimated on 
the assumption that v — 1 at the roughening transition and that the two transitions are 
separated by a thin tongue of IOF phase. 



IV. CONCLUSIONS 

We have investigated the phase transitions in a restricted solid on solid model coupled 
to an Ising model by a rejection free microcanonical Monte Carlo method. This model 
is derived from the coupled XF-Ising system by a duality transformation on the XY de- 
grees of freedom and by making what appear to be reasonable restrictions on the allowed 
height differences. The XF-Ising model is a coarse-grained generalization of the fully frus- 
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trated XY model which was originally introduced in the hope that it would make numerical 
studies of the FFXY system easier. Despite this convoluted chain of transformations and 
approximations, the RSOS-Ising model is believed to undergo continuous phase transitions 
in the same universality classes as the FFXY system. We believe that the restriction on 
the heights \Ah\ = 0, 1 is not a serious one in that it will not affect the universality class 
of the transitions provided they are continuous and that the topology of the phase diagram 
is also not qualitatively affected. We believe that the only possible difference lifting the 
RSOS restriction on the heights could have would be to allow the Ising and roughening lines 
to merge into a segment of first order transition line separating the DOF and IOR phases. 
In our numerical study of the RSOS-Ising system we identified three phases: Ising ordered 
rough (IOR), Ising ordered flat (IOF) and Ising disordered flat (DOF) which correspond 
to XY and Ising ordered, XY disorder and Ising order and complete disorder respectively. 
Using the rejection free microcanonical method, we obtain a precise phase diagram and 
demonstrate that the IOR and DOF phases are separated by a sliver of IOF phase for all 
physical values of the parameters. As the temperature in the FFXY model is raised, the 
system will follow some path in the A, C plane from the IOR phase to the DOF phase. Since 
it must pass through the intervening IOF phase, this is strong evidence that there are two 
continuous transitions in the FFXY and XF-Ising models. 

The critical exponents along the Ising transition line seem to change continuously in the 
region where the Ising and roughening transitions are close together but this may be due 
to large corrections to scaling caused by the nearby roughening transition. In that case the 
correct determination of the critical behavior will require very large systems. Since we have 
located the precise phase transition lines the true critical behavior along the Ising line may 
be determined using other MC methods appropriate for the very large systems needed. We 
leave this for future study. 
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FIGURES 

FIG. 1. Representation of bond variables by integer height variables on the dual lattice. Under 
the height variable scheme all constraints on bond variables are automatically satisfied. 

FIG. 2. (a) Surface flatness Mr of the RSOS part for lattice size L = 22. It is easily seen 
that there is a roughening transition between flat(MR ~ 1) and rough phases(MR ~ 0.5). (b) 
Magnetization M\ of the Ising part for lattice size L = 22. There is a relatively sharp transition 
region in the Ising magnetization. In the region of small A and C the Ising transition line is close 
to the roughening transition line. It should be noted that the viewing angle of (b) is 90° rotated 
from that of (a). 

FIG. 3. Specific heat Cr of the RSOS part for lattice size L = 22. The specific heat ridge 
parallel to (7-axis does not increase for L > 18 just as in the XY model, (b) Specific heat C\ of 
the Ising part for lattice size L = 22. As expected, the specific heat ridge grows as L increases. 
The location of the Ising ridge is very close to that of the RSOS ridge in small A and C region. 

FIG. 4. (a) Extrapolating the peak position of the Mr susceptibility in the bulk limit for 
several C values. Each line corresponds to C = 0.10,0.15,0.20 and 0.25 from bottom to top. 
The data points on the ^4-axis = 0) indicate extrapolated values, (b) Phase diagram of the 
RSOS-Ising model in the A, C plane. The solid and dotted lines represent the Ising and roughening 
transition, respectively. The roughening transition line is obtained by extrapolations of the Mr 
susceptibility peak positions. For the comparison we present another phase diagram in the inset, 
where the roughening transition line is the position of peak in the Mr susceptibility for L = 22. 

FIG. 5. Exponent 7/1/ of the Ising susceptibility along the Ising transition line. For A > 0.8 
7/1/ ps 1.75 as in the 2D Ising model. This value is obtained by using the best fit for the finite-size 
scaling relation for the susceptibility x ^ L" 1 ^ with L= 6, 10, 14, 18, and 22. 
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